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The .Radio Science Systems Group (RSSG) provides various support functions 
for severai Right project radio-science teams. Among these support functions are 
uplink and sequence planning , real-time operations monitoring and support, data 
validation, archiving and distribution functions , and data processing and analy- 
sis. This article describes the support functions that encompass radio-science data 
performance analysis. The primary tool used by the RSSG to fulfill this support 
function is the STBLTY program set. STBLTY is used to reconstruct observable 
frequencies and calculate model frequencies , frequency residuals , frequency stabil- 
ity in terms of Allan deviation, reconstructed phase, frequency and phase power 
spectral density, and frequency drift rates. In the case of one-way data, using an 
ultrastable oscillator (USO) as a frequency reference , the program set computes the 
spacecraft transmitted frequency and maintains a database containing the in-Right 
history of the USO measurements . The program set also produces graphical dis- 
plays. Some examples and discussion on operating the program set on Galileo and 
Ulysses data will be presented. 


I. Introduction 

The radio science investigations of the different flight projects require a diverse set of tools for ensuring 
that acquired radio science data satisfy experimenter requirements. For example, the gravitational wave 
experiments on Galileo, Ulysses, Mars Observer, and Cassini each has specific stability requirements on 
the received signal frequency and amplitude. The Radio Science Systems Group (RSSG) utilizes various 
software tools to verify that the received data satisfy the requirements and that the required configurations 
are in place. STBLTY is the primary analysis tool used by the RSSG to characterize the performance 
of the radio science instrument. The STBLTY program set has been utilized by the RSSG to perform 
measurements of ultrastable oscillator (USO) frequency for the Galileo Solar Redshift Experiment [1,2], to 
characterize the performance of the Galileo Orbiter USO [3], to verify performance of the Mars Observer 
USO, 1 and to verify performance for the gravitational wave experiments on Galileo [2] and Ulysses [4]. 
The characterization of the Galileo USO is important for occultation radio science experiments that will 
be performed during the Jupiter tour [5]. This article will focus on describing the algorithms that are 
performed by the individual component programs of the STBLTY program set and will present examples 
using STBLTY on Galileo and Ulysses radio science data. 

The version of STBLTY used for Voyager radio science was inherited by the Galileo and Ulysses 
radio science projects. This version was used primarily for one-way data occultation analysis, where 


1 D D. Morabito, “Mars Observer 93-049 USO Test Processed Through STBLTY,” JPL Interoffice Memorandum 
3394-93-077 (internal document), Jet Propulsion Laboratory, Pasadena, California, June 17, 1993. 
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a baseline frequency signature was used to remove all effects other than those due to the planetary 
atmosphere. Hence, the inherited version did not have the stringent accuracy requirements needed for 
Galileo and Ulysses radio science experiments. In order to make STBLTY suitable for Galileo and 
Ulysses cruise data, several upgrades were made to the program set after the Voyager era (Neptune 
1989). These upgrades included adding the capability of computing coherent two-way Doppler residuals, 
modeling the effects of a spinning spacecraft, upgrading the troposphere calibration model, introducing 
an ionosphere calibration algorithm using tracking system analytic calibration (TSAC) data as input, and 
estimating a more accurate trajectory model. The new trajectory model includes effects due to nutation, 
irregularities in Earth rotation (UT1-UTC), polar motion, equation of the equinoxes, zonal tides, and 
more comprehensive formulation of state vectors. 


Figure 1 is a block diagram displaying the interconnections between the component programs of the 
STBLTY program set. Input radio science data consist of closed-loop data from the DSN tracking system 
in the form of archival tracking data files (ATDFs) or open-loop data from the radio science system in 
the form of original data record (ODR) tapes. The input celestial reference set (CRS) file from the 
Navigation Team of a supported project provides the spacecraft-state vectors necessary to model the 
spacecraft trajectory from the observable data. 


INPUTS 


OUTPUTS 



Fig. 1. The STBLTY program set block diagram displaying program names, input, output, 
and intermediate tapetfile products. 

The program TRAJVECT reads in body-centered spacecraft state vectors from the CRS file, translates 
the Earth-centered vectors to the location of the observing station, and performs light-time solutions. 
TRAJVECT produces a set of downlink state vectors at the time of signal reception (in file F45), and if 
data are two-way, produces a set of uplink vectors at the time of signal transmission (in file F49). These 
data sets are passed to the program RESIDUAL. 

The closed-loop tracking system utilizes a phase-locked loop (PLL) that performs signal acquisition, 
lock-up, and detection in real time. The closed-loop data are conditioned by the Multi-Mission Navigation 
Team, which produces and delivers an ATDF to the RSSG. Among the quantities on an ATDF are Doppler 
counts, Doppler “pseudoresiduals” (residuals based on predicted frequencies used to tune the receivers), 
signal strengths (AGCs), and Doppler reference frequencies either in the form of a constant frequency 
or uplink ramps. OBSERVE reads the Doppler extractor cycle counts and Doppler extractor reference 
information from the ATDF and reconstructs these into received sky frequencies (in file F50). If data are 
two-way, OBSERVE writes the uplink ramps into the programmable ramp file (PRF). 
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The open-loop system mixes an incoming intermediate frequency (IF) signal with a signal whose 
frequency is a linear approximation of the predicted frequency, also known as the programmable oscillator 
control assembly (POCA) ramps. The baseband signal is passed through a filter whose bandpass is 
centered at the expected frequency and has a sufficiently wide bandwidth to allow for any unexpected 
signal frequency excursions. A set of analog-to-digital converters (ADCs) digitizes the received bandwidth 
and then writes the samples onto the ODR tapes. The program PLLDEC reads the ODR tapes and 
performs signal detection on the recorded samples, outputs an F36 file with the detected frequencies and 
signal strengths, and extracts the open-loop receiver tuning information (POCAs) from the ODR headers 
and writes these to an F33 file. The program OBSERVE reads the PLLDEC F36 and F33 output files 
and reconstructs the information into received sky frequencies and writes these in an F50 file that is 
passed to RESIDUAL. 

The program RESIDUAL reads the observed sky frequency F50 file and, if the data are two-way, the 
PRF file output from OBSERVE and the spacecraft trajectory vectors from the TRAJVECT F45 (one- 
way or two-way) and F49 (two-way only) output files. RESIDUAL computes model frequencies from the 
F45 and F49 trajectory vectors and a spacecraft spin model. These model sky frequencies are subtracted 
from the observed sky frequencies to produce residual frequencies that are written onto the RESIDUAL 
output F52 file. 

The F52 residual frequencies are input to the STBLTY program, which performs further corrections 
and produces Allan deviations, frequency and phase power spectral densities, and other quantities. In 
the case of one-way USO data, STBLTY computes the spacecraft transmitted frequency and writes out 
records for each pass onto a database summary file, SUMFIL, which is accessed by other programs for 
USO measurement analysis. The SUMFIL is also delivered to the radio science investigator, who performs 
an analysis for the purpose of measuring the solar gravitational redshift [1]. 


II. Program Descriptions 

A more detailed description of each of the component programs of STBLTY follows. 

A. PLLDEC — Open-Loop Data Signal Detection 

The scheme employed for detecting signal frequency and power from the digitized samples of the 
open-loop system is a digital second-order phase-locked loop program called PLLDEC. The digitized 
samples output from the ADCs in the Deep Space Communications Complex (DSCC) spectrum processing 
(DSP) subsystem are written onto the ODR tapes (see Fig. 1) that are input to PLLDEC. A theoretical 
background that supplements the following discussion is given in [6], and a detailed description of a 
similar version of PLLDEC (using noncoherent AGC requiring high signal-to-noise ratio (SNR) signals) 
is contained in an internal memorandum. 2 

PLLDEC reads the set of digitized samples, x(i),i = 1, ■ • • , TV, from the ODR tape. An initial fre- 
quency, / 0 , for the PLL is obtained from a fast Fourier transform (FFT) performed on an initial subset 
of samples. The output signal frequency, y(nT), is expressed as a linear combination (recursive function) 
of the recorded samples, rr(i), and the previous output signal-frequency estimates, y(iT). 

The implemented phase-locked loop consists of a mixer (phase detector) and loop filter providing the 
forward gain path, F(s) = (1 + as)/bs ; a low-pass filter, AT(s), for suppressing the mixer sum product 
(with negligible baseband response); and an integrator, G(s) = 1/s. The closed-loop transfer function 
(the transform of the PLL phase over the transform of the input phase) is given by 


2 A. Densmore, “A Digitally Implemented Phase- Locked Loop Detection Scheme for Analysis of Phase and Power Stability of 
a Calibration Tone,” JPL Interoffice Memorandum Voyager RSST-8S-016 (internal document), Jet Propulsion Laboratory, 
Pasadena, California, February 4, 1988. 
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L(s) 


F(s)G(s) _ 1 + as 

1 + F(s)G{s) 1 + as 4- bs 2 


The desired s-domain transfer function of the second-order PLL of the output detected frequency to the 
input phase is given by 


W(s) = 


s(l + as) 

1 + as + bs 2 


The continuous time analysis of the discrete time PLL is justified by assuming that the product of the 
PLL noise bandwidth and sample period is much less than unity, ( BpllT « 1). In order to realize the 
digital PLL implementation, the bilinear transformation from s-space to z-transform space is employed: 


It has been shown that this transformation outperforms others for the case of low sample rates [7]. 

By substituting Eq. (1) for s into G(s), the following z-domain expression for the integrator is obtained: 


G(z) — 


T 

2 



( 2 ) 


By substituting Eq. (1) for s into F(s), the corresponding z-domain function for the loop filter is obtained: 


T + 2a [1 + z~ l ((T - 2a) /(T + 2a))] 
2b 1 -z- 1 


(3) 


The relations for a and b are given by 


R + l 

4 Bpll 


R 

where R is the loop damping parameter, which is set equal to 2, the optimum value of loop performance 
specified by Jaffee and Rechtin [8] when the initial phase is unknown but uniformly distributed. A 
typical value for the one-sided loop bandwidth used for Galileo and Ulysses radio-science data processing 
is Bphh = 1 H z - 


a = 


b = 


A single-pole low-pass filter is also employed to suppress the mixer sum product and has the following 
s-domain representation: 


K(s) = 


1 

1 T (s/a) 
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Using the bilinear transformation, the implemented z-domain transfer function of the low-pass filter is 
given by 


2 + aT [1 — z - 1 ((2 - aT )/( 2 + oT))] (4) 

where a is set equal to 10 B PLL . The bilinear form of Eqs. (2), (3), and (4) are used to determine the 
appropriate coefficients that multiply the previous input samples and previous output values to yield each 
output estimate of signal phase or frequency. 

The contribution to the Allan deviation due to additive white Gaussian noise (AWGN) over the PLL 
bandwidth B P i P can be estimated from the following equation: 


<?y( T ) = 


2-KfoT 


where B PLL is the PLL bandwidth, Hz; CNR is the carrier signal-to-noise ratio of the signal; /„ is 
the observing frequency, Hz; and r is the time interval, s. By knowing the expected signal strengths 
and system noise temperatures, one can use this equation to estimate the B PLL required to obtain an 
estimate of the AWGN contribution to the Allan deviation. By reducing B P ll appropriately, one can 
reduce this contribution below other effects, thereby allowing visibility of those effects. Typical CNRs for 
Galileo and Ulysses result in AWGN contributions to the Allan deviations at 1000 s, which lie well below 
the measured values that are limited by media (for two-way Galileo and Ulysses) and USO (for Galileo 
one-way) when using 1 Hz for B PPL . 

B. TRAJVECT — Trajectory Vector Conditioning 

The program TRAJVECT reads spacecraft state vectors from the CRS file, corrects Earth-center 
vectors to the observing-station location, and performs light-time solutions to yield a set of downlink 
state vectors at the time of signal reception, and (if two-way) performs a light-time solution to yield a 
set of uplink state vectors at the time of signal transmission. The conditioned set of state vectors is then 
passed to RESIDUAL. For discussions on how the NAV orbit determination solutions are performed, 
from which the CRS files are derived, the reader is referred to [9] for Galileo and [10] for Ulysses. 

1. Input CRS File State Vectors. For cruise processing, where the gravitational effects of the 
planets are assumed negligible, it suffices that the CRS file contains EME 1950 (Earth mean equator and 
equinox of B1950.0) state vectors of the spacecraft relative to the Sun and the Earth, as follows: 

X ® -8 /c(f/t) = Earth center-to-spacecraft position vector 
v® a / c (tk) = Earth center-to-spacecraft velocity vector 
X 0“ s / c (tfc) = Sun-to-spacecraft position vector 
V 0-8 /c(tfc) = Sun-to-spacecraft velocity vector 

where fc is the index of the vectors which are evenly spaced in time, T s apart. The time tag, t k , is 
referenced at the spacecraft in ephemeris time in seconds past 1950.0. 

2. Conversion to Sun Center. The state vectors from the CRS file are referred to Sun center using 
the following translations: 
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X®-®(tfc) = 

x®- s / c (f fe )-x 0 - 8 / c (t fc ) 

(5) 

x s/ c -0(^) = 

- x 0 - 8 / c (t fc ) 

(6) 

v®-°(tfc) = 

V0_S/C(tfc) (Jei) ~ v °" S/C( * fc) 

(7) 

v s / c - 0 (tk) = 

- v°- a/c (tk) 

(8) 


where the (+/rel) notation means that the velocities add relativistically. 

3. Calculation of Observing Station Vectors. TRAJVECT translates the Earth-Sun vectors 
[Eqs. (5) and (7)] to the location of the observing station. The station-Earth vectors are computed in 
the EME 1950 frame and are added to the Earth-Sun vectors [Eqs. (5) and (7)] to produce the required 
station-to-Sun state vectors. Among the corrections applied to the station vectors are precession, nutation, 
UT1-UTC, polar motion, equation of the equinoxes, and solid Earth tides. 

A table of station locations in cylindrical coordinates (referenced to the mean pole, equator, and prime 
meridian of 1900.05) are read to obtain the observing station coordinates: radius off the spin axis, r s 
(km); height above the equatorial plane, 2 (km); and longitude in degrees east of the Greenwich meridian, 
A. Thus, the body-fixed station vector in rectangular coordinates is 


x d " 


r s cos A 

r 3 sin A 

z 


Several notations are used to denote time: the corresponding Julian date of t*, 

tl D = — ■ + 2,433,282.5 

k 86,400 

the number of Julian centuries of 36,525 mean solar days that have elapsed since 1950.0, T k , 


T k 


tk 

3,155,760,000 


and the difference of ephemeris time (ET) and universal time (UT1), A t k - 

To translate the station vectors to the frame of the current epoch, several translations and rotations 
are applied. Polar motion matrix rotations are applied, and the vectors are then rotated to the correct 
position in the frame of the current epoch. The rotation matrices that perform these corrections use values 
of universal time and polar motion obtained from the monthly publications of the International Earth 
Rotation Service (IERS). These values are interpolated to the appropriate time tag, t k , of the vectors, 
yielding the Earth orientation values of At UT iR (UT1R-UTC) and polar motion X p and Y p . The A turiR 
is further corrected by adding the effects of short-period zonal tides (monthly and fortnightly), to yield 
Atari • The zonal tides are computed using the formulation of Yoder et al. [11]. 

The angle of the Greenwich meridian at the current epoch, corrected for variation in Earth rotation, 
ac(tk), is known as Greenwich mean sidereal time, or the Greenwich hour angle of the mean equinox of 
date, and is given by 
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a G(tk) — C*o(£fc) + a\{tk) + a 2 (<fc) + a 3 (<fc) 


(9) 


where a 0 {t k ) = -(1.002737909294 + T fc 0.589 x 10- ]O )((At fc )/240) 4- 100.0755426042 + 1^0.0003870833 is 
the fraction of a turn (in degrees) of the Earth over T k Julian centuries, ai(t k ) = MOD[(t k / 240), 360] is 
the fraction of a turn (in degrees) due to daily rotation, a 2 (t k ) = MOD[36,000T fc ,360] is the fraction of 
a turn (in degrees) due to the revolution of the Earth in its orbit, and a 3 (t k ) = 0.76931208337*; is the 
fraction of a turn (in degrees) over T k Julian centuries, which accounts for the slowing down of the Earth 
due to the exertion of tidal forces. 

The nutation model used by the program set is the International Astronomical Union (IAU) 1980, 
or Wahr nutation model, which was adopted by the IAU in 1984 and uses a series developed by Wahr 
based on the nutation model of P. K. Seidelmann [12]. The expected accuracy of this model (1 milliarc) 
[13], translates to about a 1-mHz error at 8.4 GHz, but is expected to be constant over a tracking pass. 
The angles of nutation estimated from this model are Ac, the nutation in obliquity; A T , the nutation in 
longitude; and e m , the mean obliquity of the ecliptic. The derivatives are computed from the differences of 
the nutation angles evaluated at two slightly different times divided by the time difference. The obliquity 
is added to the mean obliquity to yield the true obliquity of the ecliptic, e. 

The Greenwich mean sidereal time in Eq. (9), a G {t), is converted to Greenwich apparent sidereal time, 
a GA{t) (the Greenwich hour angle of the true equinox of date), by applying the correction of the equation 
of the equinoxes: 


a GA(t k ) = a G (t k ) + A'l' cose 


The space-fixed station vector projected in the celestial frame at t k (the position of date) is computed 
from the body-fixed station vector by applying a series of rotations that corrects for polar motion (R t and 
R 2 ) and Earth rotation (R 3 ). This “true of-date” space-fixed station vector is corrected for nutation (N), 
which performs the rotation from the true of-date equator and equinox to the “mean of-date” equator and 
equinox. Finally, the station vector is rotated from the mean of-date equator to the mean equator of the 
reference epoch (EME 1950) using the precession rotation matrix, P. See the Appendix for definitions of 
P and N. The rotation matrices Ri,R 2 , and R 3 are given by 


R3(ao»(4)) = 


cos a GA (t k ) -sin a GA (t k ) 0 
sin a GA (t k ) cos a G A (t k ) 0 
0 0 1 


R 2 (X P ) = 


COS X p 
0 

sin X p 


0 

1 

0 


- sin X p 
0 

cosXp 


Rift) 


1 0 0 

0 cos Yp sin Y p 
0 - sin Y p cos Y p 


The station position vector in EME 1950 coordinates is computed as follows: 


x d “-®(t fc ) = PN[R 3 R 2 R 1 x dss + Ax tide .] 


( 10 ) 
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where Ri, R 2 , and R 3 are the “Earth platform” rotation matrices given above, N is the nutation rotation 
matrix, P is the precession rotation matrix, and A Xtides(tk) designates the perturbations due to solid 
Earth tides which, along with velocity perturbations, Av tides, are induced by the Sun and Moon. 3 

The station velocity vector is given by 4 ' 5 ’ 6 

v dss ~®(t k ) = [PN + PN][R 3 R 2 Rix dss + A xudes{t k )} + PN[R 3 R 2 Rix dss + Av tirfes (t fc )] (11) 


where the N and P matrices (see the Appendix) denote the derivatives of the nutation and precession 
matrices, respectively, A Vtidesi.t'k') the station velocity perturbations due to tides, and R 3 is given by 


R-a^CA^fc)) = 


-sin a GA {tk ) 
cos a G A{tk) 

0 


- cos a GA {Tk) 0 

— sin a G A(tk) 0 

0 0 


4. Light-Time Solutions. The light-time solution is performed as follows. First, an initial estimate 
of the one-way light time is made using Eqs. (5) and (6): 


To(tfc) = 


lx 0 ~ e (tfc) - x»/ c -°(t fc )| 
c 


where the “0” subscript denotes the initial estimate. The ith estimate of the time tag of downlink 
reception a one-way light time later is made as follows: 


tk,i ~ h + n(tk,i-l) (tk,o = tk) 


( 12 ) 


The three sets of vectors with time tags, and tj+i, lying closest to the estimated received time, 

tk,u are quadratically interpolated to yield estimates of the state vectors at t k ,i- The station-Earth 
vectors [Eqs. (10) and (11)] are evaluated at t fcil and added to the Earth-Sun vectors (Eqs. (5) and (7) 
interpolated at £/e,t) as follows: 


x dss-0 


v dss-0 


= X dSS_ ®(ffc,i) +X® _e (ffe,i) 

(13) 

(tfc.i) = v dss -®(f fcli ) v®-©(t fcli ) 

(14) 


The correction to the light-time of the next iteration is 


3 T. D. Moyer, “Correction to Earth Fixed Station Coordinates Due to Solid Earth Tides, Ocean Loading, and Pole Tides 
and Calculation of Periodic Terms of UT1,” JPL Engineering Memorandum 314-505 (internal document), Jet Propulsion 
Laboratory, Pasadena, California, July 4, 1991. 

4 N. A. Mottinger and T. D. Moyer, “A Close Encounter With a Doppler Observable,” JPL Interoffice Memorandum 
314.5-1025 (internal document), Jet Propulsion Laboratory, Pasadena, California, June 18, 1984. 

5 T. D. Moyer, “Proposed Changes to ODP Transformation Between Body-Fixed and Space-Fixed Coordinates for the 
Planets and the Sun,” JPL Engineering Memorandum 314-271 (internal document), Jet Propulsion Laboratory, Pasadena, 
California, June 16, 1982. 

6 T. D. Moyer, “Time Derivative and Partial Derivatives of the Body-Fixed to Space-Fixed Coordinate Transformation for 
the Sun, Planets and Planetary Satellites,” JPL Engineering Memorandum 314-379 (internal document), Jet Propulsion 
Laboratory, Pasadena, California, August 22, 1985. 
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Ar i+1 (t M ) = 


|x dss - Q (^, i )-x s / c -Qfa, < )| 


7~i(tk,i— l) 


(15) 


c 


and the new light-time is 


r i+i(4,i) = Ti(t kt i-i) 4- Ar i+1 (* fc>i ) 


The time tag for the next iteration is refined as [Eq. (12), substituting i 4 1 for i] 




The correction [Eq. (15)] is compared against a preset tolerance that should be tight enough for 
accurate vector determination given the stability requirements of the experiment being processed and the 
time spacing of the input CRS vectors. The current value used for Galileo and Ulysses cruise data is 
10~ 4 s, assuming vectors separated by 30 s. Once the criterion is met, the iterative solution is completed 
and returns the estimated vectors at If the criterion is not met, then additional iterations are 

performed until it is satisfied. In the case of the uplink vectors, the procedure is the same except that a 
minus sign replaces the plus sign in the appropriate equations. 

The time tag of the final iteration of the light-time solution for the downlink output vectors is redefined 
for notational convenience as 




The M state vectors from the CRS file are processed this way such that a set of M heliocentric station 
and spacecraft state vectors referenced at times fj, k = 1, ■ • ■ , M is produced in the EME 1950 reference 
frame. These vectors are passed to RESIDUAL via the downlink vector file F45 and uplink vector file 
F49. For downlink, 

r 0 (*:) 
v£r 0 (*z) 

v£' 0 «> 

For uplink, 

<?-°(*?) = station position 
v up S ~ G (t*) = station velocity 
Xup C °(tj T) = spacecraft position 
Vup C G (t*) = spacecraft velocity 

where the subscript k is the index of the kth vector in the F45 downlink vector file, and the subscript / 
is the index of the Zth vector in the F49 uplink vector file. 


= station position 
= station velocity 
= spacecraft position 
= spacecraft velocity 
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C. OBSERVE — Sky Frequency Reconstruction 

1. Closed-Loop System. The program OBSERVE reads in “Doppler” cycle counts from the closed- 
loop system, (f>(tj), along with information needed to compute the Doppler extractor reference frequency, 
fref(tj )■ The cycle counts consist of an integer number of counts output from the Doppler counter, plus 
a fractional term output from the Doppler resolver. The Doppler reference frequency may either be a 
constant synthesizer frequency (SIM), normally used as the reference during one-way data acquisition, 
or be computed from uplink ramps used as the reference during two-way data acquisition. The ramps 
used to tune the uplink frequency consist of a series of programmable frequencies and ramps that drive 
a digitally controlled oscillator (DCO), which is multiplied up and input to the exciter that drives the 
transmitter. The ramps are written to the PRF file, which is input to RESIDUAL, where they are used 
to reconstruct the uplink frequencies for two-way data. 

For closed-loop data read off an ATDF for a standard Deep Space Station (DSS), OBSERVE computes 
the received biased Doppler frequencies from the accumulated cycle counts as follows: 


r \ 0(+l-l) ~ 0(*j-l) 

Jdop\ l j) “ , . 

tj + l Cj_i 


In the case of a modulo reset, where <f>(tj) is less than 1 ), 


fdop{tj) 






-u 


The Doppler frequency is converted to a 2.3-GHz sky frequency by 


= 96— frefitj) -B (f dop (t 3 ) - 10 6 ) 


where 10 6 is the 1-MHz bias, B is the sign of the bias, f Te f is the Doppler extractor reference frequency- 
interpolated to the time tag of the measured Doppler, and 240/221 is the spacecraft transponder ratio. 
For 8.4 GHz, 


/*(*;) = 96— / re/ + ) -B (f dop ( tj) - 10 6 ) 

In the case of the 34-m high-efficiency (HEF) station, the 8.4-GHz frequency is given by 7 * 

f aky (tj) = ^^[32(4. 68125/ r e/ (tj) - 81.4125 10 6 ) + 6500 10 6 ] -B (/ dop (t/) - 10 6 ) 

2. Open-Loop System. For the case of open-loop data, the frequencies detected from the recorded 
samples by PLLDEC, / rec , in the F36 file are converted to sky frequencies, f 9 or /*, using the tuning 
frequencies, / poca , interpolated from the POCA offsets and ramps in the F33 file at the f rec time tag. 
For the radio science IF-VF converter assembly (RIV) receivers, 


7 T. D. Moyer, “Change to the ODP and ODE for Processing X-Band Uplink Data,” JPL Engineering Memorandum 314-430 

(internal document), Jet Propulsion Laboratory, Pasadena, California, October 15, 1987. 
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f s (tj) = 3 (fpocaitj) + (™ x 10 6 ^ + (1950 x 10 6 ) - f samp + f rec (tj) 

f x (tj) = 11 (fpocaitj) - (10 x 10 6 ) + (8050 X 10 6 ) - 3 f samp + f rec (tj) 
where f 3am p * s the sampling frequency of the recorded samples. For the multimission receivers (MMRs) 


S s (tj) = 4 8f poca (tj) + (300 


10 6 ) f 

1 w samp 


frec{tj ) 




n 

T 


48 fpocaitj) + (300 X 10 6 ) 


If 

A J samp 


“f freci^j) 


The observed sky frequencies from OBSERVE are input to RESIDUAL, where model and residual fre- 
quencies are calculated. 

D. RESIDUAL — Residual Frequency Calculation 

The program RESIDUAL reads the sky frequencies from the OBSERVE output F50 file and the 
state vectors from the TRAJVECT output F45 and F49 files to compute model frequencies based on 
the spacecraft trajectory, and produces frequency residuals for either 2.3 GHz, 8.4 GHz, or differenced 
S-3/11X data. RESIDUAL applies a polarization-dependent correction for a spinning spacecraft 8 and 
can estimate and remove a sinusoidal signature induced by an off-axis antenna on a spinning spacecraft 
(e.g., low-gain antenna (LGA)-2 on Galileo). 

1. Spin Correction. For the constant spin correction, specific default values are built into the code 
to account for different spacecraft configurations. This correction assumes the 2.3-GHz one-way default 
values for f spin shown in Table 1 (note that Galileo has a 3-rpm spin rate and Ulysses has a 5-rpm 
spin rate). This correction could also be obtained from attitude-control telemetry data and the defaults 
therefore overridden. For Galileo passes involving LGA-2, the spin rate could also be solved for from a 
three-parameter sinusoid fit of the trajectory corrected residuals. 

Table 1. One-way default values for f spxn . 


Spacecraft 

Antenna 

Spin mode 

fS,\-way H 
* spin ’ 

Galileo 

LGA1, HGA 

Dual 

0.052 


LGA2 

Dual 

-0.052 


LGA1, HGA 

All 

0.048 


LGA2 

All 

-0.048 

Ulysses 



-0.0833 

Any other 



0.00 


RESIDUAL applies the appropriate spin correction to the sky frequency depending upon the spacecraft 
spin configuration and frequency band. If the data are 2.3-GHz one-way, f spin is applied as given in 
Table 1. For the case of two-way data with a 2.3-GHz downlink and a 2.3-GHz uplink, 


8 P. Priest and J. Breidenthal, “Galileo Spin Effects,” JPL Interoffice Memorandum 3393-90-186 (internal document), Jet 
Propulsion Laboratory, Pasadena, California, November 30, 1990. 
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fS/S 

J spin 



rSil-way 
J spin 


For the case of two-way data with an 8.4-GHz downlink and a 2.3-GHz uplink, 


f s/x _ 

J spin 


240 11 , 

— + 1 

221 3 


rS, 1 — way 
f spin 


For the case of dual-frequency S-X residuals, the spin correction applied is given by 

fS-X _ ° rS,l-way 

J spin 11 *^ 

2. One-Way Data Processing. The state vectors from the TRAJVECT F45 file are interpolated 
to the time tag of the observable frequency {tf ss at the observing station and t x at the spacecraft) and 
are used to estimate the downlink Doppler correction factors: 


d i 


(c - |ys/c-0(^/ c )|)( c+ |y S /c-6(^/ c )|) 


/ c \ 


(16) 


v dBS-© ( t ds*) x ds 9 -0 ( t d»B) _ x s/c -0( t V c ) 

c |x dss -©(«f ss ) - x s / c -©(t, s/c )| 

v"/e-©(tj/ c ) X dss - 0 (tf s ) - x s/c ~ e (f- /C ) 
C |x dss -©(tf s ) - x s / c - 0 (ti /c )| 


(17) 


(18) 


d 4 



v ds3-0( t dBB)|)[ c _)_ 


|v d38 ~0(t ds3 )|] 


(19) 


The gravitational redshift correction factors that account for the frequency shifts due to the masses of 
the Sun and the Earth, evaluated at f t , are given by 


G Mp, 


|j x dss-©| | x »/c-©|J 


dm — 


G Mr 


[|x d8S -®| |x s / c “®|J 


The radiated 2.3-GHz frequency of the spacecraft is computed from the received frequency of the first 
data point as follows: 


f s/cii i) 


{did^/ di,di) + d® + d® 


+ /, 


spin 
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where the Doppler and redshift terms are evaluated at t A . The estimated sky frequencies (evaluated at t t ) 
are computed from this transmitted spacecraft frequency at the first data point (evaluated at 0 using 
the above correction factors (evaluated at f t ) as follows: 


f S (U)=fs/c(t l) 


d\d,2 


+ d Q + 


For the one-way case, the residuals are defined such that the first data point is anchored to zero. Given 
the set of observable frequencies output from the OBSERVE F50 file, the estimated frequency for 

the first data point is, by definition, 




The residual frequencies are computed by differencing the model from the observable as follows: 


A f s (ti) = f 3 (t i )-f*(t l ) 

The residual frequencies and the spacecraft transmitted frequency are passed to STBLTY in the F52 file 
for further processing. 

3. Two-Way Data Processing. This section discusses the computations of two-way Doppler 
frequencies and residuals. The downlink Doppler corrections [Eqs. (16)— (19)] are evaluated at the time 
tag of the received signal frequency, tf ss , as described in the previous section. In addition, the uplink 
vectors from the F49 file are interpolated to the time of uplink signal transmission, t“ p , and the Doppler 
correction terms for the uplink signal leg are similarly computed as follows: 


Ui(t.) = u 

/ (c - |v d88 “0(^ p )|)(c + |v dss -©(^ p )|) 

( 20 ) 

c 2 

u 2 {t t ) = 1 

V»/c-©(^/c) x »/c- 0 (^/c) _ X d»«-©^UP) 

C | X a /°-©(t’ /c ) _ x d,.- 0(t « P) | 

( 21 ) 

U 3 (ti) = 1 

V daa -©(t“ p ) x a / c -©(«’ /c ) - x daa -0(t“ p ) 
C |x B /c-0( ( a A) _x daa -©(t“ P )| 

( 22 ) 

u 4 (*i) = y 

1 (c - |v a /«=-©(t* /c )|)(c + |v a /'-©(t* /c )|) 

(23) 

c 2 


The radiated frequency at the time of signal transmission, is obtained from the uplink ramps in the 
PRF file, interpolated at and then multiplied up to sky frequency, The estimated 2 . 3 -GHz 

received frequency is computed as follows: 


f 3 (U) 



240 U 1 U 2 

221 U3U4 


f s/s 


d$d4 
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The 2.3-GHz residual frequency is computed from the observed and estimated received frequencies at 
time ti as follows: 


a nu) = f 3 (t i )-r(t l ) 


Likewise, for 8.4 GHz, 

d 1^2 

d-^d. i 

A f x (t i ) = f x (t i )-f x (t i ) 


f x {U) = 


ft (fP) ±1 240 mu 2 

U ; 3 221 u 3 u 4 


rs/x 
* spin 


4. Correction for Off- Axis Antenna on a Spinning Spacecraft. A sinusoidal signature due to 
an off-axis antenna on a spinning spacecraft can be removed by fitting a three- parameter model to the 
bias-corrected residuals: 


A f 3 (t l ) = C l sin (C 2 ti + C 3 ) 


where U is the time tag of the data point. The spacecraft transmitted frequency is corrected so that it 
is referenced to the phase center of the spacecraft spin. This is done by removing the value at the first 
data point: 


fs/c(t 1 ) = fs/c(t 1 ) - Ci sm(C 3 ) 


The residuals are adjusted accordingly: 


A f 3 c {U) = A f(ti) + Ci sin {C 2 ti + C 3 ) - Ci sin(C 2 ti + C 3 ) 

The spacecraft transmitted frequency, / s / c , is passed to the program STBLTY in the F52 file along with 
other header information, observed frequencies, and residuals. 

E. STBLTY — Media Correction, Display Plots, and Output File 

STBLTY is the final program that applies additional corrections to the data, writes out records to a 
database file for USO and gravitational redshift analysis, and produces plots. At this point, the residual 
frequencies A f s {t t ) and/or A f x (U) for L data points i = 1, 2, • • • , L are read in from the RESIDUAL F52 
output file. If both bands are available, STBLTY then combines the observable frequencies to produce 
the differential data type: 


A r(U) = - T f x (U ) (24) 

For one-way data, the spacecraft transmitted frequency referenced to the first data point, f s / c {t i)> 
passed in the header. 

1. Ionosphere/ Charged-Particle Calibration. The charged particle effect can be “removed” 
from the radio science data in one of two ways: (1) When dual- frequency data are available, a differential 
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correction computed from the differential data type, Eq. (24), is applied to the individual frequency bands 
or (2) a TSAC-supplied polynomial describing the ionosphere path delay signature can be converted into 
frequency and removed from the frequency residuals. The TSAC correction removes ionospheric effects 
for either dual- or single-band downlink passes. 

For the case where two simultaneous downlink frequency channels exist, the charged particle can be 
removed using the differential frequency type f s - 3/ x /ll, given in Eq. (24). Given that the frequency 
residuals to be corrected are A f s (t t ) for 2.3 GHz and A f x (t t ) for 8.4 GHz, the corrected residuals are 


A/ C 5 (ti) = A/ a (*i) - ~Af sx {ti) 

for a one-way 2. 3- GHz downlink; 

33 

Af-(t i ) = Af x (t i )-—Af sx (U) 


for a one-way 8.4-GHz downlink; 


A/c fat) = A - 


'240 

221 


121 121" 
112 + U2 


A /“(tO 


for a two-way 2.3-GHz uplink/2. 3-GHz downlink; and 


A /*(«,) 


A - 


'240 

221 


121 33 

112 + Il2 


A /“(«*) 


for a two-way 2.3-GHz uplink/8.4-GHz downlink, where the subscript c denotes the corrected value of 
the residual frequency. 

The TSAC group produces and delivers ionosphere calibration files to the RSSG, where they are 
archived and sometimes delivered to the radio science experimenters. The files contain polynomial coeffi- 
cients describing the mapped ionosphere path delay in meters (referenced to a 2295-MHz carrier) for each 
DSCC and for each time period that overlaps a tracking pass for a given spacecraft. TSAC determines 
these coefficients from observations of Earth-orbiting satellites at each complex and maps to the line of 
sight of the spacecraft. 

The TSAC coefficients and time tags are converted to frequency corrections in Hz (referenced to 
2.3 GHz) by STBLTY. Given the start time, (£*,), and end time, (£/), of the interval for which the 
polynomial is defined, and the polynomial coefficients (a u i = 0, ■ • ■ , 5), the ionosphere path delays to the 
spacecraft at regularly spaced intervals of time tj are computed as follows: 


5 

a Pion(tj) - y a, j* 

t=0 


where 


Xj 


= 2 


lj_ ^6 

if — 


- 1 
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The path delays computed at time offsets ±T / 2 about t, are differenced and scaled to yield the one-way 
2.3-GHz ionosphere frequency correction at tj-. 


,, , f a [A p im (tj + (T[ 2)) - Apian (tj - (r/2))] 

/ ion \rj ) 

Figure 2 is an example of a typical path delay profile derived from a TSAC polynomial. Figure 3 is the 
corresponding ionosphere frequency profile at 2.3 GHz. 

The set of time tags, tj, and corrections, / ion (tj), that overlaps the tracking pass is used to perform 
the calibration only if dual-frequency charged-particle calibration was not previously performed. Given 
the frequency residual, A f 3 {t t ) or A/ I (t i ), and the above ionosphere correction interpolated to time 
t l ,f- on (t l ), the corrected residuals are then computed as follows: For a one-way 2.3-GHz downlink, 


A/*(ti) = A/ 3 (tj) - fi on {ti) 


for a one-way 8.4-GHz downlink, 


A/ C *(ti) = A r(t t ) - ^f° on (U) 


for a two-way 2.3-GHz uplink/2. 3-GHz downlink, 


240 

A/c (ti) - A f(U) - fioniU) - ~ *r) 


and for a two-way 2.3-GHz uplink/8. 4-GHz downlink, 


3 11 240 

A/*(t0 = A r(U) - —fioniU) - J 221 f!on(tr - U) 


where t r is the round-trip light time. 

For the c^e of differential f s - 3/ x /ll two-way data, the contribution of the common 2.3-GHz uplink 
will be zero after differencing the two downlink channels. The remaining charged-particle contribution 
in the A f sx data type (SX) will be due to the downlink only. Therefore, this data type is corrected as 
follows: 


A /r(ti) = A/“(ti) - (25) 

2. Troposphere Calibration. The troposphere calibration is performed after the ionosphere cali- 
bration so that the remaining time-dependent signature of the residuals over elevation angle more closely 
follows the troposphere ray-path function. The troposphere can be removed using either a model with a 
specified zenith troposphere path delay (default 2.1-m) or the zenith path delay can be fit and removed 
from the data if there is sufficient arc over elevation angle. The ray-path mapping function used is the 
CfA (Center for Astrophysics) mapping function [14]: 


136 


FREQUENCY, Hz 



Fig. 2. An example of a typical ionosphere path delay profile derived from 
a TSAC polynomial. 



Fig. 3. An example of a typical frequency profile over a pass derived from 
the TSAC delay profile of Fig. 2. 
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1 


sin e(ti) + (bi/ (tan e(i) + (6 2 / (sin e(U) + 63))) 

where e(U) is the elevation angle and the coefficients b k are functions of meteorological parameters, given 

by 


b\ = 0.0002723[1 + (2.642 x 10 -4 p 0 ) - (6.400 x 10“ 4 e 0 ) + (1.337 x 10 2 T D ) - (8.550 x 10 2 a) 

- (2.456 x 10- 2 /i 2 )] 

b 2 = 0.0004703(1 + (2.832 x 10 _5 p 0 ) + (6.799 x 10~ 4 e 0 ) + (7.563 x 10" 3 T 0 ) - (7.390 x 10 _2 a) 

- (2.961 x 10 _2 /i 2 )] 

63 = -0.0090 

Reasonable values of surface pressure (p 0 = 900 mbar), partial water vapor pressure (e 0 = 5 mbar), 
surface temperature (T 0 = 292 K), tropopause altitude (h 2 = 12.2 km), and dry model parameter 
(a = 5.0) are used in the computation of these coefficients. This model is in agreement with predictions 
made by Lanyi [15] to the 1-cm level down to a 5-deg elevation angle. 

The one-way troposphere correction to the observed 2.3-GHz frequency residual is given by 


f tropic) 


/ a (t,)Ap z 

cT 


' / T\ 

n f T\ 

R ( t{ 4* ~ 

— R [ti — - 

L V 2 ) 

V 2/. 


where / s (tj) is the downlink sky frequency, A p z is the zenith troposphere path delay, c is the velocity 
of light, and Tj 2 is a chosen time offset about f t . For the cases of two-way 2.3-GHz uplink/2. 3-GHz 
downlink and 2.3-GHz uplink/8. 4-GHz downlink, the corrections are 


ftrop(ti) 


-/ 5 (t,)Ap; 

cT 


R ( ti + 


*H)] + s[ r K-‘0 


- R I ti - — - t, 


ftropiU) 


-/ I (t i )Ap z 

cT 


R [ ti 4* 




T 

2 


where t r is the round-trip light time. The one-way frequency residuals are adjusted accordingly: 

A fHU) - - fl rov {ti) -f / t a rop (ti) 


A f x c {U) = A f x (ti) - f* rop (U) + f? rop (ti) 


effectively removing the signature of the troposphere relative to the first data point, because f 3 / c (t 1) 
was evaluated here. A bias correction evaluated from the average of the corrected residuals is applied 
to the transmitted spacecraft frequency to refer it to the center of mass of the residuals. The spacecraft 
transmitted frequency at 2.3 GHz is corrected as follows: 


fc,s/c(tmid) — fs/c(t 0 + ftropft l) + A 


where the last term is the average of the corrected residuals over the data span and is the time tag 
of the midpoint of the data acquisition interval. 
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The two-way frequency residuals are corrected for troposphere as follows: 


a/e(*i) = a nu) - r trop { tl ) 
= Af x (ti) - ftropHi) 


III. STBLTY Performance Analysis Results and Discussion 

The following sections will discuss examples using STBLTY on Galileo and Ulysses flight data. 

A. Galileo Results 

Examples of running STBLTY on Galileo radio science data will be presented in the following two 
sections. 

1. One-Way Data. The results of routine processing of Galileo USO data using STBLTY will be 
briefly discussed in this section. For further details on this analysis of closed-loop one-way USO data 
acquired during the first 2 years of cruise (1989-1991), the reader is referred to [3]. 

Due to the undeployed high-gain antenna (HGA), only 2.3-GHz data from the LGA antennas were 
available. The passes were conducted weekly on the average, were typically 2 hours in duration, and were 
centered about meridian crossing where the troposphere effect was minimal. For troposphere calibration, 
it was sufficient to remove a model (versus fit and remove) using a default zenith path delay. The error 
incurred by not performing an ionosphere calibration is expected to be at the mHz level or below for a 
2-hour pass centered about meridian crossing at 2.3 GHz. 

The Allan deviation results were consistent with the values expected at the appropriate time intervals, 
given the known signal levels and noise mechanisms. At the short time intervals, the noise in the signal 
was dominated by thermal noise due to the low signal levels that were due to transmitting through 
the LGA. At the long time intervals, the noise in the signal was dominated by the USO itself. The 
transmitted USO-referenced spacecraft frequency followed a behavior that was consistent with known 
aging mechanisms of the USO, having an rms scatter of 17 mHz about a fitted aging model. 

2. Two-Way Data. The testing of different calibration schemes using STBLTY on Galileo two-way 
2.3-GHz data will be discussed in this section. STBLTY was run on three passes: (1) a DSS-63 tracking 
pass on day 91-123, (2) a DSS-63 pass on day 93-081, and (3) a DSS-14 pass on day 93-082. The latter 
two passes are from the 1993 Gravitational Wave Experiment. 

Table 2 displays the processing results for these three passes using different calibration schemes. The 
first column denotes the year and day number of the pass; the next column is the DSN station identifica- 
tion; the next column denotes the ionosphere calibration employed (see the key for an explanation of the 
codes); the next column denotes the troposphere calibration employed (see the key for an explanation of 
the codes); the next column is the zenith path delay and uncertainty in the troposphere calibration (if 
the code in the TRP column is “F,” this value is from the least-square fit); the next column is the slope 
from a linear fit of the residuals after calibrations have been performed; and the last four columns are 
the Allan deviations of the postcalibration residuals at 1, 10, 100, and 1000 s. 

Four different calibration schemes were employed for each pass. The first run (ION:T, TRP:F) involved 
calibrating the ionosphere utilizing TSAC polynomial coefficients and calibrating the troposphere by 
fitting the model over the elevation angle for a zenith path delay and bias. The second run (ION:N, 
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TRP:F) involved no ionosphere calibration and calibrating the troposphere as described in the first 
scheme. The third run (ION:T, TRP:M) involved calibrating the ionosphere using TSAC polynomial 
coefficients and calibrating the troposphere by removing the model using the default zenith path delay of 
2.1 m (no fitting performed). The last run (ION:N, TRP:N) did not employ any media calibration. 

The zenith path delay for the TSAC ionosphere/troposphere fit run (T,F) for passes 91-123 and 93-081 
are in reasonable agreement with the expected 2.1-m zenith path delay. The zenith path delay fit from the 
data when no ionosphere is removed (N,F) is consistently lower than those of the TSAC calibrated runs 
(T,F) and (T,M). This suggests that the TSAC calibration does contribute to removing nontropospheric 
trends from the data. 

The linear slope of the residuals tends to increase in magnitude as fewer calibrations are performed. 
The magnitude of the linear slope tends to be smaller for the runs involving both troposphere fit and 
TSAC ionosphere calibration (T,F) and larger for the run involving no calibrations (N,N). 

The Allan deviations at the time intervals of 1, 10 and 100 s do not differ significantly between the 
different schemes over each given pass. The Allan deviation at 1000 s for the three runs involving a tropo- 
sphere calibration (T,F), (N,F), and (T,N1) are in reasonable agreement. The 1000-s Allan deviations are 
consistent with values predicted due to solar plasma at the known solar elongation angles using the model 
of Armstrong, Woo, and Estabrook [16]. However, the Allan deviations at 1000 s are significantly larger 
for the no-media calibration case (N,N), clearly showing that the removal of the systematic troposphere 
signature with elevation angle is important for obtaining meaningful results at this time interval. 

For Galileo at 2.3 GHz, trajectory errors of 0.018-mm/s in velocity translate to Allan deviations at 
1000 s of about 6 x 10“ 14 . This lies below the 10“ 13 level expected due to plasma (for two-way data) and 
the instability of the USO (for one-way data). 


Table 2. Galileo media calibration test results. 


Pass, 

yr-day 

DSS 

ION 

TRP 

Zenith path 
delay, 
m 

Linear 
slope, 
10 -7 Hz/s 

1(T 14 

v„(10), 
10" 14 

ITy(lOO), 

icr 14 

<Tj,{1000), 

icr 14 

91-123 

63 

T 

F 

1.96 ± 0.03 

0.12 

± 

0.29 

— 

357 

27.7 

24.8 

91-123 

63 

N 

F 

0.98 ± 0.03 

0.16 

± 

0.29 

— 

357 

27.6 

27.0 

91-123 

63 

T 

M 

2.1 

0.84 

± 

0.29 

— 

357 

27.7 

25.0 

91-123 

63 

N 

N 

— 

-5.01 

± 

0.29 

— 

357 

27.7 

38.2 

93-081 

63 

T 

F 

2.17 ± 0.03 

0.77 

± 

0.27 

1947 

206 

25.4 

17.9 

93-081 

63 

N 

F 

1.70 ± 0.04 

1.55 

± 

0.27 

1947 

206 

25.5 

23.3 

93-081 

63 

T 

M 

2.1 

0.58 

± 

0.27 

1947 

206 

25.4 

18.4 

93-081 

63 

N 

N 

— 

-3.55 

± 

0.27 

1947 

206 

26.1 

50.5 

93-082 

14 

T 

F 

1.11 ± 0.04 

0.39 

± 

0.32 

1808 

194 

26.9 

14.5 

93-082 

14 

N 

F 

0.80 ± 0.05 

0.75 

± 

0.31 

1808 

195 

26.9 

17.9 

93-082 

14 

T 

M 

2.1 

2.73 

± 

0.32 

1808 

195 

26.9 

13.4 

93-082 

14 

N 

N 

— 

-1.14 

± 

0.32 

1808 

195 

27.1 

29.1 


T = TSAC ionosphere calibration. 

D = Dual- frequency (f s - 3/ x /ll) charged-particle calibration. 

N = No calibration performed. 

F = Troposphere calibration (model removed using fit zenith path delay). 

M = Troposphere calibration (model removed using default zenith path delay). 
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B. Ulysses Results 

For Ulysses, two sets of two-way data were analyzed: (1) closed-loop data from a dual-band pass 
from the first opposition 91-004 (DSS 14) and (2) open-loop data from a dual-band pass from the second 
opposition 92-064 (DSS 14). Table 3 displays the results for these data sets and for each of the applicable 
data types within each pass. For the single-band data (2.3 and 8.4 GHz), the following runs were 
performed: (1) TSAC ionosphere and troposphere fit (T,F), (2) dual-frequency calibration of charged 
particles, and troposphere fit (D,F), (3) no ionosphere/charged-particle and troposphere fit (N,F), and 
(4) TSAC ionosphere and troposphere model removed using the default zenith path delay (T,M). The 
no-calibration case (N,N) is not presented in this table since the (N,N) Galileo runs (Section III. A) clearly 
illustrate the degradation when not applying a media correction. The description of Table 3 is the same 
as that of Table 2 in the proceeding section except that the data type code (S, X, or SX) follows the 
year-day number in the first column. 


Figure 4 is a plot of the 8.4-GHz residuals for pass 91-004 for the fully calibrated (T,F) run. Figure 5 
is a plot of the averaged residuals where every 60 points have been averaged to allow long period trends 
to be easily examined. Figure 6 is a plot of the reconstructed phase for this pass, and Fig. 7 is a plot of 
the corresponding Allan deviation. Figure 8 is a plot of the averaged residuals prior to the removal of the 
troposphere, thus illustrating the signature of an unmodeled troposphere. Figure 9 displays the elevation 
profile over the pass and the troposphere model that was fit from the data of Fig. 8 and then removed 
from the residuals, resulting in the plots displayed in Figs. 4-7. Figure 10 displays the TSAC ionosphere 
correction at 2.3 GHz for this pass. 


Table 3. Ulysses media calibration test results. 


Pass, 

yr-day 

DSS 

ION 

TRP 

Zenith path 
delay, 
m 

Linear 
slope, 
10 -7 Hz/s 

10~ 14 

<r w (10), 

10-14 

a v (100), 

icr 14 

<Ty(l000), 

IQ-14 

91-004 S 

14 

T 

F 

1.50 ± 0.06 

-0.001 ± 0.71 

1080 

116 

15.4 

4.8 

91-004 S 

14 

D 

F 

1.90 ± 0.08 

0.08 ± 1.71 

2624 

276 

32.5 

2.4 

91-004 S 

14 

N 

F 

0.44 ± 0.06 

0.09 ± 0.71 

1081 

116 

15.4 

6.2 

91-004 S 

14 

T 

M 

2.1 

1.35 ± 0.71 

1081 

116 

15.4 

6.7 

91-004 X 

14 

T 

F 

1.69 ± 0.03 

0.29 ± 1.65 

621 

66 

9.3 

1.9 

91-004 X 

14 

D 

F 

1.40 ± 0.04 

0.65 ± 1.83 

756 

80 

12.8 

3.7 

91-004 X 

14 

N 

F 

1.05 ± 0.04 

0.42 ± 1.65 

621 

66 

9.3 

2.8 

91-004 X 

14 

T 

M 

2.1 

3.51 ± 1.65 

621 

66 

9.3 

3.6 

91-004 SX 

14 

T 

N 

— 

0.66 ± 0.67 

1107 

116 

13.3 

3.3 

91-004 SX 

14 

N 

N 

— 

1.51 ± 0.67 

1107 

116 

13.4 

5.2 

92-064 S 

14 

T 

F 

1.54 ± 0.05 

0.19 ± 0.12 

335 

46 

16.7 

6.7 

92-064 S 

14 

D 

F 

1.65 ± 0.04 

-0.004 ± 0.12 

384 

45 

16.4 

7.4 

92-064 S 

14 

N 

F 

1.17 ± 0.04 

0.13 ± 0.12 

335 

46 

16.7 

5.9 

92-064 S 

14 

T 

M 

2.1 

1.26 ± 0.12 

335 

46 

16.7 

6.6 

92-064 X 

14 

T 

F 

1.43 ± 0.03 

0.12 ± 0.25 

147 

30 

13.6 

5.6 

92-064 X 

14 

D 

F 

1.46 ± 0.04 

0.15 ± 0.27 

187 

31 

13.8 

6.1 

92-064 X 

14 

N 

F 

1.38 ± 0.04 

0.23 ± 0.25 

147 

30 

13.6 

5.7 

92-064 X 

14 

T 

M 

2.1 

4.63 ± 0.24 

147 

30 

13.6 

6.2 

92-064 SX 

14 

T 

N 

— 

0.15 ± 0.09 

303 

32 
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Fig. 9. The (a) elevation angle profile and (b) troposphere correction fit from the frequency 
residuals of Fig. 8. This signature was then removed from the data, resulting in the 
calibration illustrated in Figs. 4 through 7. 
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Fig. 10. Ionosphere frequency correction at 2.3 GHz, derived from a TSAC polynomial 
for the DSS-14 Ulysses pass of January 4, 1991. 


The fitted zenith path delays for case (T,F) appear reasonable but are somewhat low compared to 
the expected 2.1-m values. Unmodeled trends, which blend with the signature of the troposphere, are 
assumed to exist in the data and thus influence the fit. 

The linear slope from a straight line least-squares fit of the postcalibrated residuals tends to be of 
insignificant or of marginal magnitude when a troposphere fit is performed. When a model troposphere 
is removed, the fitted linear slopes of the residuals are significant. This suggests that significant unknown 
trends may be present in the data set that are absorbed by the troposphere fit or that the troposphere 
model using the default zenith path delay is not sufficient. 

The Allan deviations at 1, 10 and 100 s do not change significantly between the different runs for a 
given band, except for the dual-frequency charged-particle calibrated runs (D in the ION column), where 
significant system noise at these time intervals is introduced into the data. At 1000 s, dual-frequency 
charged-particle calibration (D,F) results in the lowest Allan deviations at 2.3 GHz for the 91-004 pass 
and no significant changes at 2.3 GHz for pass 92-064. At 8.4 GHz, the TSAC calibrated run for pass 
91-004 produced the lowest Allan deviation, 1.9 x 10' 14 . For the 92-064 DSS-14 pass, no significant 
improvement was noted in the 8.4-GHz 1000-s Allan deviation between the different runs. The 91-004 
Allan deviations agree with independently measured values [4]. 

For the f 9 - 3/ x /ll data type (SX), the following runs were performed: (1) TSAC correction 

[Eq. (25)](T) and no troposphere (N) and (2) no ionosphere calibration (N) and no troposphere cali- 
bration (N). Since the difference frequency eliminates all nondispersive effects, troposphere calibration 
is not performed. For 91-004 SX and 92-064 SX, the slopes of the calibrated runs were of insignificant 
or marginal magnitudes, suggesting the effectiveness of using TSAC polynomials to remove long period 
trends from the SX data type and illustrating the dominance of the time variability of the ionosphere 
over that of the solar plasma during these observations near solar opposition. The Allan deviations 
do not change significantly between runs for the SX data type, except that the 1000-s value for pass 
91-004 SX shows an improvement from 5.2 x 10“ 14 to 3.3 x 10' 14 when applying the TSAC calibration. 
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Figure 11 displays the averaged residuals for this pass without applying TSAC (N,N). Here, all nondis- 
persive effects should cancel out, and the resulting signature is presumably due to dispersive sources. The 
corresponding TSAC S-X frequency calibration for this pass is shown in Fig. 12. The averaged residuals 
after applying the TSAC calibration (T,N) are shown in Fig. 13. The resulting signature is considerably 
flatter than that of the non-TSAC calibration run (N,N) of Fig. 11, but a long period trend of order 
0.5-mHz over this pass is apparent upon closer inspection. This is consistent with the observation of 
Bertotti et al. [4] that the error of the TSAC calibration is not negligible and can be as large as 0.5 mHz. 
The long period trends on the order of 0.5 mHz that remain in the data after applying media calibrations 
are consistent with known uncertainties in the calibration schemes. The zenith path delays for the TSAC 
ionosphere can be in error by as much as 75 cm for day passes and 15 cm for night passes. The troposphere 
can be in error by 5 cm. 9 * An overall 15-cm zenith path delay error (for these nighttime solar opposition 
passes) can map into variations at the 1-mHz level over a pass, consistent with what is observed in these 
data. 

For Ulysses at 8.4 GHz, trajectory errors of 0.012 mm/s in velocity at 1000 s translate to Allan 
deviations at about 4 x 10” 14 , which is comparable with the observed values. These are also comparable 
with expected noise due to media. After removal of the ionosphere and troposphere calibration models, 
significant random fluctuations in the data remain. These fluctuations have different characteristics with 
different frequencies, solar elongations, local weather conditions, and local ionospheric variability. These 
are usually masked by system noise at the lower time intervals, whose magnitudes depend upon received 
signal levels and received noise power of a given spacecraft. Expected levels of these fluctuations at 1000 s 
for Ulysses at 8.4 GHz are about 3 x 10” 14 for the troposphere, an upper bound of 3 x 10” 14 for plasma 
noise, and 2.5 x 10” 14 for the system noise of the closed-loop receiver [4]. 

Figure 14 displays the 8.4-GHz Allan deviation curves for a set of different calibration schemes for the 
91-004 DSS-14 pass. For the case of two-way data over a sufficient elevation angle arc, the best estimate 
of system stability is obtained by applying a TSAC ionosphere correction and by fitting and removing a 
troposphere model. 


IV. Conclusion 

A description of the STBLTY program set has been presented with details on model implementation 
and residual analysis. Examples of the use of this program set on Galileo and Ulysses data were also 
presented. The frequency of a spacecraft signal received on the ground is the observable data type. All 
known effects are removed from this observable to produce residuals that can be inspected to evaluate 
performance. Among the errors in the residuals are those due to white noise, media, and trajectory. For 
a sufficiently high SNR, random media effects appear to be the limiting error source at the 1000-s time 
scale for both Galileo and Ulysses coherent data sets. 

Future upgrades to STBLTY will focus on modeling the gravity fields of planets in order to remove 
these effects from the data during flybys or orbital operations. The first use for this capability will involve 
the Galileo Jupiter orbital operations when Galileo is expected to tour the Jovian system starting in late 
1995. Other future upgrades include modifying the code to process 32.0-GHz downlink data such as that 
expected for the Cassini radio science experiments, to process coherent downlink data that involve 8.4- 
or 32.0-GHz (Ka-band) uplink signals, converting the code to manipulate vectors in the J2000 system, 
and incorporating a planetary atmospheric model for occupation data analysis. The program set was 
recently modified to process three-way data (implemented in a developmental version). 


9 T. McElrath, personal communication, Navigation Systems Section, Jet Propulsion Laboratory, Pasadena, California, 

1994. 
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Fig. 11. Averaged differenced frequency residuals (dual-band) for the DSS-14 
Ulysses pass of January 4, 1991. 



Fig. 12. Ionosphere frequency correction for dual-band derived from a TSAC 
polynomial for the DSS-14 Ulysses pass of January 4, 1991. 
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Fig. 14. Allan deviation curves for a set of different calibration schemes for 
the DSS-14 Ulysses pass of January 4, 1991. 
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Appendix 

Formulation of Precession and Nutation Matrices 


The following definitions for system rotation matrices and their derivatives apply when defining the 
precession and nutation rotation matrices and their derivatives: 
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I. Precession Rotation Formulation 

The precession matrix, P, is the product of the following three system rotations using the definitions 
given above: 


P “ (2 ” a 1950^ - <5l95o) ( — A 195 o) 2 

where the precession rotation angles in terms of Julian centuries since 1950, T, are given by 

A 1950 = [89.9999986565 - 0.64027801T - (3.042075 x 10~ 4 )T 2 - (5.0837 x 10" 6 )T 3 1 — 

' v ^ J 180 
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The derivative of the precession matrix is as follows: 
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[-0.64027801 - (2 x 3.042075) x 10 _4 T - (3 x 5.0837) x 10_6t2 ]^ 
[-0.5567500297 + (2 x 1.185607) x 10~ 4 T + (3 x 1.16119) x 10 _5 T 2 ]^ 
[-0.6402780091 - (2 x 8.39481) x 10~ 5 T - (3 x 0.50003) x 10_5t2 ]^ 

II. Nutation Rotation Formulation 

The nutation rotation matrix is the product of the following three rotations using the matrix rotation 
definitions given above: 

N = (-e m )x(A*) a (e m + Ae), 

where the nutation angles are the IAU 1980 model values defined in the text. The derivative of the 
nutation rotation matrix is thus given as 

N = - « («m + £) (A tf) 2 (€ m + Ae) x + AV{-e m ) x (-AV + ^ 

x (tm + Ac)* + e(— e m ) x (A ^') 2 e m — Ae + — ^ 


dt ‘ = 
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j7 a 1950 - 
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